/***************************************************************************
 *   Scientific Construct (sC++)						                   *
 *   The Scientific Computational Toolkit				                   *
 *   Discrete Geometry Module                                              *
 *   Copyright (C) 2006 by Shekhar S. Chandra           		           *
 *   Shekhar.Chandra@sci.monash.edu.au                                     *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; Version 2 of the License.               *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#ifndef DGRADONTRANFORM_H
#define DGRADONTRANFORM_H

#include <blitz/tinymat.h>
#include <blitz/array.h>

using namespace blitz;

/*! \class DGRadonTransform
    \brief Class designed to contain and encapsulate Discrete Radon Transforms.
        The transforms can be operated on Blitz Array or Matrices.
    @author Shekhar S. Chandra
*/
template <typename type, int size, int rank>
class DGRadonTransform
{
public:
    // -------------------------------------------------------
	/*! \name Constructors */
	//@{
    /*! \fn DGRadonTransform::DGRadonTransform(void)
	   	\brief Standard Constructor
	   */
    DGRadonTransform(void);
    /*! \fn DGRadonTransform::~DGRadonTransform(void)
	   	\brief Standard Destructor
	   */
    virtual ~DGRadonTransform(void);
    //@}

    // -------------------------------------------------------
	/*! \name Radon Transforms */
	//@{
    /*! \fn DGRadonTransform::FiniteRT(Array<type,rank> &field)
        \brief The Finite Radon Transform of Matus & Flusser (1993). It is
            assumed that the field size is prime. Overloaded for Blitz Array.
        */
    Array<type,rank> FiniteRT(Array<type,rank> &field);
	/*! \fn DGRadonTransform::FiniteRT(Array<type,rank> &field, Array<type,rank> & bins)
        \brief The Finite Radon Transform of Matus & Flusser (1993). It is
            assumed that the field size is prime. Overloaded for Blitz Array.
        */
    void FiniteRT(Array<type,rank> &field, Array<type,rank> &bins);
    void FiniteRT_modulo(Array<type,rank> &field, Array<type,rank> &bins, int modulo);
    /*! \fn DGRadonTransform::iFiniteRT(Array<type,rank> &field)
        \brief The Inverse Finite Radon Transform of Matus & Flusser (1993). It is
            assumed that the field size is prime. Overloaded for Blitz Array.
        */
    Array<type,rank> iFiniteRT(Array<type,rank> &bins);
	/*! \fn DGRadonTransform::iFiniteRT_Unscaled(Array<type,rank> &field)
        \brief The Inverse Finite Radon Transform of Matus & Flusser (1993). It is
            assumed that the field size is prime. Overloaded for Blitz Array.
        */
    Array<type,rank> iFiniteRT_Unscaled(Array<type,rank> &bins);
	/*! \fn DGRadonTransform::iFiniteRT(Array<type,rank> &bins, Array<type,rank> &result);
        \brief The Inverse Finite Radon Transform of Matus & Flusser (1993) done to recover only a subset of the image. It is
            assumed that the field size is prime. Unlike other functions, size of result is assumed to be set to size of subset.
            Overloaded for Blitz Array.
        */
    void iFiniteRT(Array<type,rank> &bins, Array<type,rank> &result);
    void iFiniteRT_modulo(Array<type,rank> &bins, Array<type,rank> &result, Array<type,1> &inverses, int &modulo, bool norm = true);
	/*! \fn DGRadonTransform::iFiniteRT_Unscaled(Array<type,rank> &bins, Array<type,rank> &result);
        \brief The Inverse Finite Radon Transform of Matus & Flusser (1993). It is
            assumed that the field size is prime. Overloaded for Blitz Array.
        */
    void iFiniteRT_Unscaled(Array<type,rank> &bins, Array<type,rank> &result);
    /*! \fn DGRadonTransform::iFiniteRT_Subset(Array<type,rank> &bins, Array<type,rank> &result);
        \brief The Inverse Finite Radon Transform of Matus & Flusser (1993) done to recover only a subset of the image. It is
            assumed that the field size is prime. Unlike other functions, size of result is assumed to be set to size of subset.
            Here the x and y coordinates represents the position of the top left-hand corner of the image.
            Overloaded for Blitz Array.
        */
    void iFiniteRT_Subset(Array<type,rank> &bins, Array<type,rank> &result, const int xcoord, const int ycoord);
    void iFiniteRT_Subset_Unscaled(Array<type,rank> &bins, Array<type,rank> &result, const int xcoord, const int ycoord);
    /*! \fn DGRadonTransform::FiniteRT(TinyMatrix<type,size,size> &field)
        \brief NOT OPERATIONAL. The Finite Radon Transform of Matus & Flusser (1993). It is
            assumed that the field size is prime. Overloaded for Blitz Matrix. Returns matrix of bins. NOT OPERATIONAL
		@todo Fix ISum Problem in Inverse.
        */
    TinyMatrix<type,size+1,size> FiniteRT(TinyMatrix<type,size,size> &field);
    /*! \fn DGRadonTransform::iFiniteRT(TinyMatrix<type,size+1,size> &bins)
        \brief NOT OPERATIONAL. The Inverse Finite Radon Transform of Matus & Flusser (1993). It is
            assumed that the field size is prime. Overloaded for Blitz Matrix. Returns reconstructed matrix. NOT OPERATIONAL
        */
    TinyMatrix<type,size,size> iFiniteRT(TinyMatrix<type,size+1,size> &bins);
    //@}

    // -------------------------------------------------------
	/*! \name I/O */
	//@{
	/*! \fn DGRadonTransform::init()
        \brief Initializes all internal variables.
        */
	inline void init()
	{   Isum = initializer;		}
	/*! \fn DGRadonTransform::setIsum(type sum)
        \brief Sets the \f$I_{sum}\f$ calculated in RT's (i.e. the total image sum).
        */
	inline void setIsum(type sum)
	{   Isum = sum;				}
	/*! \fn DGRadonTransform::getIsum()
        \brief Returns the \f$I_{sum}\f$ calculated in RT's (i.e. the total image sum).
        */
	inline type getIsum()
	{   return Isum;            }
	/*! \fn DGRadonTransform::ISumFromBins(Array<type,rank> &bins)
        \brief Returns the \f$I_{sum}\f$ calculated in RT's (i.e. the total image sum).
        */
	void ISumFromBins(Array<type,rank> &bins);
	/*! \fn DGRadonTransform::getInitializer()
        \brief Sets the initializer value for the RT's.
        MUST BE SET FOR CORRECT RESULTS WHEN USING TYPES OTHER THAN NUMBER TYPES (eg. string).
        */
	inline void setInitializer(type value)
	{   initializer = value;            }
	/*! \fn DGRadonTransform::getInitializer()
        \brief Returns the initializer value for the RT's.
        MUST BE SET FOR CORRECT RESULTS WHEN USING TYPES OTHER THAN NUMBER TYPES (eg. string).
        */
	inline type getInitializer()
	{   return initializer;            }
    //@}

protected:
    type Isum; /*!< \brief The sum of all the Field values */
    type initializer; /*!< The initialial value to which the RT's are to be conducted */
};

template <typename type, int size, int rank>
DGRadonTransform<type,size,rank>::DGRadonTransform()
{

}

template <typename type, int size, int rank>
DGRadonTransform<type,size,rank>::~DGRadonTransform()
{

}

template <typename type, int size, int rank>
Array<type,rank> DGRadonTransform<type,size,rank>::FiniteRT(Array<type,rank> &field)
{
    const int p = field.cols(), translate_x = 1;
    int translate_y, y, k;
    Array<type,rank> bins(p+1,p);
    Isum = initializer;
    bins = initializer;

    //Get Projections where j denotes projection number
    for(int j = 0; j < p; j ++)
    {
        y = 0;
        for(int m = 0; m < p; m++)
		{
			y += j;
			k = (m + translate_x)%p;
			for(translate_y = 0; translate_y < p; translate_y ++)
				bins(j,translate_y) += field(k, (translate_y + y)%p);
		}
    }
    //Get the last projection of rows
    for (int j = 0; j < p; j ++)
        for(k = 0; k < p; k ++)
        {
            bins(p,j) += field(j,k);
            Isum += field(j,k);
        }

    return bins;
}

template <typename type, int size, int rank>
void DGRadonTransform<type,size,rank>::FiniteRT(Array<type,rank> &field, Array<type,rank> &bins)
{
	const int p = field.cols(), translate_x = 1;
    int translate_y, y, k;
    bins.resizeAndPreserve(p+1,p);
    Isum = initializer;
    bins = initializer;

    //Get Projections where j denotes projection number
    for(int j = 0; j < p; j ++)
    {
        y = 0;
        for(int m = 0; m < p; m++)
		{
			y += j;
			k = (m + translate_x)%p;
			for(translate_y = 0; translate_y < p; translate_y ++)
				bins(j,translate_y) += field(k, (translate_y + y)%p);
		}
    }
    //Get the last projection of rows
    for (int j = 0; j < p; j ++)
        for(k = 0; k < p; k ++)
        {
            bins(p,j) += field(j,k);
            Isum += field(j,k);
        }
}

template <typename type, int size, int rank>
void DGRadonTransform<type,size,rank>::FiniteRT_modulo(Array<type,rank> &field, Array<type,rank> &bins, int modulo)
{
	const int p = field.cols(), translate_x = 1;
    int translate_y, y, k;
    bins.resizeAndPreserve(p+1,p);
    Isum = initializer;
    bins = initializer;

    //Get Projections where j denotes projection number
    for(int j = 0; j < p; j ++)
    {
        y = 0;
        for(int m = 0; m < p; m++)
		{
			y += j;
			k = (m + translate_x)%p;
			for(translate_y = 0; translate_y < p; translate_y ++)
			{
				bins(j,translate_y) += field(k, (translate_y + y)%p);
				bins(j,translate_y) %= modulo;
			}
		}
    }
    //Get the last projection of rows
    for (int j = 0; j < p; j ++)
        for(k = 0; k < p; k ++)
        {
            bins(p,j) += field(j,k);
            bins(p,j) %= modulo;
            Isum += field(j,k);
        }
}

template <typename type, int size, int rank>
Array<type,rank> DGRadonTransform<type,size,rank>::iFiniteRT(Array<type,rank> &bins)
{
    const int p = bins.cols(), translate_x = 1;
    int translate_y, y, k;
    Array<type,rank> result(p,p);
    result = initializer;

	if(Isum == 0)
		ISumFromBins(bins);

    //Get Projections where j denotes projection number
    for(int j = 0; j < p; j ++)
    {
        y = 0;
        for(int m = 0; m < p; m++)
		{
			y += p-j;
			k = (m + translate_x)%p;
			for(translate_y = 0; translate_y < p; translate_y ++)
				result(j,translate_y) += bins(k, (translate_y + y)%p);
		}
    }

	for(int j = 0; j < p; j ++)
		for(k = 0; k < p; k++)
		{
			result(j,k) += bins(p,j);
			result(j,k) = static_cast<type>( (result(j,k) - Isum)/p );
		}

    return result;
}

template <typename type, int size, int rank>
Array<type,rank> DGRadonTransform<type,size,rank>::iFiniteRT_Unscaled(Array<type,rank> &bins)
{
	const int p = bins.cols(), translate_x = 1;
    int translate_y, y, k;
    Array<type,rank> result(p,p);
    result = initializer;

	if(Isum == 0)
		ISumFromBins(bins);

    //Get Projections where j denotes projection number
    for(int j = 0; j < p; j ++)
    {
        y = 0;
        for(int m = 0; m < p; m++)
		{
			y += p-j;
			k = (m + translate_x)%p;
			for(translate_y = 0; translate_y < p; translate_y ++)
				result(j,translate_y) += bins(k, (translate_y + y)%p);
		}
    }

	for(int j = 0; j < p; j ++)
		for(k = 0; k < p; k++)
		{
			result(j,k) += bins(p,j);
			result(j,k) = static_cast<type>( (result(j,k) - Isum) );
		}

    return result;
}

template <typename type, int size, int rank>
void DGRadonTransform<type,size,rank>::iFiniteRT(Array<type,rank> &bins, Array<type,rank> &result)
{
	const int p = bins.cols(), translate_x = 1;
    int translate_y, y, k;
    result.resizeAndPreserve(p,p);
    result = initializer;

	if(Isum == 0)
		ISumFromBins(bins);

    //Get Projections where j denotes projection number
    for(int j = 0; j < p; j ++)
    {
        y = 0;
        for(int m = 0; m < p; m++)
		{
			y += p-j;
			k = (m + translate_x)%p;
			for(translate_y = 0; translate_y < p; translate_y ++)
				result(j,translate_y) += bins(k, (translate_y + y)%p);
		}
    }

	for(int j = 0; j < p; j ++)
		for(k = 0; k < p; k++)
		{
			result(j,k) += bins(p,j);
			result(j,k) = static_cast<type>( (result(j,k) - Isum)/p );
		}
}

template <typename type, int size, int rank>
void DGRadonTransform<type,size,rank>::iFiniteRT_modulo(Array<type,rank> &bins, Array<type,rank> &result, Array<type,1> &inverses, int &modulo, bool norm)
{
	const int p = bins.cols(), translate_x = 1;
    int translate_y, y, k;
    result.resizeAndPreserve(p,p);
    result = initializer;

	if(Isum == 0)
		ISumFromBins(bins);
    Isum %= modulo;

    //Get Projections where j denotes projection number
    for(int j = 0; j < p; j ++)
    {
        y = 0;
        for(int m = 0; m < p; m++)
		{
			y += p-j;
			k = (m + translate_x)%p;
			for(translate_y = 0; translate_y < p; translate_y ++)
			{
                result(j,translate_y) += bins(k, (translate_y + y)%p);
                result(j,translate_y) %= modulo;
			}
		}
    }

	for(int j = 0; j < p; j ++)
		for(k = 0; k < p; k++)
		{
			result(j,k) += bins(p,j);
			result(j,k) %= modulo;
			result(j,k) = (result(j,k) - Isum + modulo)%modulo;
			if(norm)
                result(j,k) = inverses(result(j,k)); //Divide by p
		}
}

template <typename type, int size, int rank>
void DGRadonTransform<type,size,rank>::iFiniteRT_Unscaled(Array<type,rank> &bins, Array<type,rank> &result)
{
	const int p = bins.cols(), translate_x = 1;
    int translate_y, y, k;
    result.resizeAndPreserve(p,p);
    result = initializer;

	if(Isum == 0)
		ISumFromBins(bins);

    //Get Projections where j denotes projection number
    for(int j = 0; j < p; j ++)
    {
        y = 0;
        for(int m = 0; m < p; m++)
		{
			y += p-j;
			k = (m + translate_x)%p;
			for(translate_y = 0; translate_y < p; translate_y ++)
				result(j,translate_y) += bins(k, (translate_y + y)%p);
		}
    }

	for(int j = 0; j < p; j ++)
		for(k = 0; k < p; k++)
		{
			result(j,k) += bins(p,j);
			result(j,k) = static_cast<type>( (result(j,k) - Isum) );
		}
}

template <typename type, int size, int rank>
void DGRadonTransform<type,size,rank>::iFiniteRT_Subset(Array<type,rank> &bins, Array<type,rank> &result, const int xcoord, const int ycoord)
{
    const int p = bins.cols(), translate_x = 1;
    int translate_y, y, k;
    result = initializer;

	if(Isum == 0)
		ISumFromBins(bins);

    //!\todo Bounds checking here is not accounted for. Segfault when xcoord+result.rows() > bins.rows().

    //Get Projections where j denotes projection number
    for(int j = xcoord; j < result.rows()+xcoord; j ++)
    {
        y = 0;
        for(int m = 0; m < p; m++)
		{
			y += p-j;
			k = (m + translate_x)%p;
			for(translate_y = ycoord; translate_y < result.cols()+ycoord; translate_y ++)
				result(j-xcoord,translate_y-ycoord) += bins(k, (translate_y + y)%p);
		}
    }

	for(int j = xcoord; j < result.rows()+xcoord; j ++)
		for(k = 0; k < result.cols(); k++)
		{
			result(j-xcoord,k) += bins(p,j);
			result(j-xcoord,k) = static_cast<type>( (result(j-xcoord,k) - Isum)/p );
		}
}

template <typename type, int size, int rank>
void DGRadonTransform<type,size,rank>::iFiniteRT_Subset_Unscaled(Array<type,rank> &bins, Array<type,rank> &result, const int xcoord, const int ycoord)
{
    const int p = bins.cols(), translate_x = 1;
    int translate_y, y, k;
    result = initializer;

	if(Isum == 0)
		ISumFromBins(bins);

    //!\todo Bounds checking here is not accounted for. Segfault when xcoord+result.rows() > bins.rows().

    //Get Projections where j denotes projection number
    for(int j = xcoord; j < result.rows()+xcoord; j ++)
    {
        y = 0;
        for(int m = 0; m < p; m++)
		{
			y += p-j;
			k = (m + translate_x)%p;
			for(translate_y = ycoord; translate_y < result.cols()+ycoord; translate_y ++)
				result(j-xcoord,translate_y-ycoord) += bins(k, (translate_y + y)%p);
		}
    }

	for(int j = xcoord; j < result.rows()+xcoord; j ++)
		for(k = 0; k < result.cols(); k++)
		{
			result(j-xcoord,k) += bins(p,j);
			result(j-xcoord,k) = static_cast<type>( (result(j-xcoord,k) - Isum) );
		}
}

template <typename type, int size, int rank>
TinyMatrix<type,size+1,size> DGRadonTransform<type,size,rank>::FiniteRT(TinyMatrix<type,size,size> &field)
{
    const int p = size, translate_x = 1;
    int translate_y, y, k;
    TinyMatrix<type,size+1,size> bins;
    Isum = initializer;
    bins = initializer;

    //Get Projections where j denotes projection number
    for(int j = 0; j < p; j ++)
    {
        y = 0;
        for(int m = 0; m < p; m++)
		{
			y += j;
			k = (m + translate_x)%p;
			for(translate_y = 0; translate_y < p; translate_y ++)
				bins(j,translate_y) += field(k, (translate_y + y)%p);
		}
    }
    //Get the last projection of rows
    for (int j = 0; j < p; j ++)
        for(k = 0; k < p; k ++)
        {
            bins(p,j) += field(j,k);
            Isum += field(j,k);
        }

    return bins;
}

template <typename type, int size, int rank>
TinyMatrix<type,size,size> DGRadonTransform<type,size,rank>::iFiniteRT(TinyMatrix<type,size+1,size> &bins)
{
    const int p = size, translate_x = 1;
    int translate_y, y, k;
    TinyMatrix<type,size,size> result;
    result = initializer;

    //Get Projections where j denotes projection number
    for(int j = 0; j < p; j ++)
    {
        y = 0;
        for(int m = 0; m < p; m++)
		{
			y += p-j;
			k = (m + translate_x)%p;
			for(translate_y = 0; translate_y < p; translate_y ++)
				result(j,translate_y) += bins(k, (translate_y + y)%p);
		}
    }

	for(int j = 0; j < p; j ++)
		for(k = 0; k < p; k++)
		{
			result(j,k) += bins(p,j);
			result(j,k) = static_cast<type>( (result(j,k) - Isum)/p );
		}

    return result;
}

template <typename type, int size, int rank>
void DGRadonTransform<type,size,rank>::ISumFromBins(Array<type,rank> &bins)
{
    for(int k = 0; k < bins.cols(); k ++)
        Isum += bins(0,k);
}

#endif // DGRADONTRANFORM_H
